function [varargout] = PSC_Scaled_meanvar(icur, isamp, m)

global CONTROL
global TRACE

if(nargin == 0) % get from current selection
    sf = getmainselection;
    EPSC = CONTROL(sf).EPSC2;
    icur = EPSC.icur;
    isamp = EPSC.isamp;
    m = EPSC.nevents;
    icrec = EPSC.icrec;
    icstart = EPSC.icstart;
    icwin = EPSC.icwin;
    bl = EPSC.ibl;
end;
fprintf(1, 'Array size: [%d %d], nevents: %d\n', size(icur,1), size(icur,2), m);

hmv = newfigure('Mini_analysis2', 'Mean-Var analysis');
% use the "scaled mini variance" method of Traynelis et al., 1993.
% 1. align each event along the steepest part of the rising phase (shifts
% points) 
% >>> find peak
% >>> find point of max dv/dt in smoothed waveform (S-G smoothing, heavy)
% calculate shift.
% 1b. scale peak of mean waveform to peak at same time in aligned wv
% subtract (for variance analysis)
% 2. sigma^2(I)  = i * I - I^2/ N + sigmab^2
% where sigmab is the variance of the background noise
% i is the apparent unitay conductance I is the current, and N is the
% estimated # of channels.
% 2b: the mean and variance is calculated by binning the data
%
span = 5; % moving average smoothing
window = ones(span, 1)/span;
ism = convn(icur(1:m, :)', window, 'same')'; %Smooth the data
ism = icur(1:m, :);
dism = diff(ism,1,2);
[imin, maxdI] = min(dism, [], 2); % find peak rising phase
% generate alignment
icmin = min(maxdI);
icmax = max(maxdI); % what is the minimum delay?
% for each trace, shift it to the left to match the minimum latency event
iev_mean2 = mean(icur, 1);
iev = zeros(size(icur));
k = 0;
for i = 1:m
    ishift = maxdI(i)-icmin;
    st = icstart(i)-ishift;
    iev(i,:) = TRACE(icrec(i), st:st+icwin(i)-1)-bl(i);
    [mbvar bvar(i)] = mean_var(icur(i, 1:10));
end;
% span = 5; % moving average smoothing
% window = ones(span, 1)/span;
% ism = convn(iev(1:m, :)', window, 'same')'; %Smooth the data
for i = 1:size(iev, 2)
    iev_mean(i) = mean_var(iev(:, i));
    ism_mean(i) = mean_var(ism(:, i));
end;
[ism_pk jsm_pk] = min(ism_mean); % get the position of the min of the average;
subplot(4,1,1)
plot(1:size(iev, 2), ism);
hold on
plot(1:size(iev, 2), ism_mean, 'r-', 'LineWidth', 2.0);

% now calculate the scaled mean and variance for each point in time:
npts = size(iev, 2);
ievave = zeros(npts, 1);
ievvar = zeros(npts, 1);
ievn = zeros(npts, 1);
idiff = zeros(size(iev));
nsamp = zeros(npts, 1);
subplot(4,1,2);
% ensemble variance, scaled
for i = 1:m % through each detected event
    ipk = min(ism(i,:)); % get peak from smoothed version of this trace
    idiff(i,:) = iev(i,:)-(ipk*iev_mean/min(iev_mean)); % difference from scaled mean
end;
sumsq = idiff .^2; % get sum squared differences
ievvar = sum(sumsq)/(m-1); % sum and calculate variance

% plot sca;ed variance as well. 
subplot(4, 1, 2);
plot((1:length(ievvar))+jsm_pk, ievvar);

iimin = min(iev_mean); % min current...
nbin = 20;
tmean = zeros(nbin, 1); 
tvar = zeros(nbin, 1);
binsize = iimin/nbin; % bin size for binned analysis.

for i = 1:nbin
    binbot = (i-1)*binsize;
    bintop = i*binsize;
    bindata = [];
    binsumsq = 0;
    [t0, it0] = findnear(iev_mean(jsm_pk-1:end), bintop);
    [t1, it1] = findnear(iev_mean(jsm_pk-1:end), binbot);
    eventsinbin = [it0:1:it1]+jsm_pk-2;
    tvar(i) = mean(ievvar(eventsinbin));
    tmean(i) = mean(iev_mean(eventsinbin));
end;

[tmean; tvar]

subplot(2,1,2);
plot(-tmean, tvar-mean(bvar), 'ko', 'markersize', 2.0); % , ievvar, 'ko', 'markersize', 2.0);
% 2. sigma^2(I)  = i * I - I^2/ N + sigmab^2
% where sigmab is the variance of the background noise
% i is the apparent unitay conductance I is the current, and N is the
% estimated # of channels.
knan = find(~isnan(tmean) & ~isnan(tvar));
[P] = polyfit(-tmean(knan), tvar(knan)-mean(bvar), 2);
A = P(1);
B = P(2);
sb2 = P(3);
g_sc = B;
N = -1/A;
%Po = ievave(2)*(-P(1)/P(2));
fprintf(1, 'N: %6.1f  gamma = %7.3f   sigmab: %6.3f  (bvar: %6.3f)\n', ...
    N, g_sc, sb2, mean(bvar));
imppl = 0:1:ceil(max(abs(tmean)));
ppl = polyval(P, imppl);
hold on
plot(imppl, ppl, 'r-');
